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ABSTRACT 



A method is developed, using the quadratic form of canonical un- 
coupled state variables as a Lyapunov function to determine the switch- 
ing logic for a quasi-optimum control of a non-linear dynamic system 
Analytio arguments are presented to support the method, which is capable 
of being extended to any order system and any number of controls, sub- 
jeot to oertain practical limitations noted in the analysis . 

A computational scheme for determining the canonical state variables 
and control functions is presented and a topological interpretation is 
made of the choice of variables used for the control logic calculation.. 

The controller based on the method described is applied to a non- 
linear second-order system. Phase trajectories for both the uncontrolled 
and controlled systems are obtained by means of a digital computer simula- 
tion. Various aspects of the theoretical limitations of the method pre- 
sented are investigated and the experimental results are analyzed with 
respect to the theoretical predictions. 

#The system is postulated to be described by a system of differential 
equations of known form. 



ii 



TABLE OF CONTENTS 



Section Title Page 

1 0 Introduction 1 

2. Theory 5 

3 e Application of Theory to Design of a Single 

Controller 16 

4# Simulation of Control of a Non-Linear Second 

Order System 20 

5 e Conclusions 36 

Appendix I - Flow Charts and FORTRAN Language 
Program for Digital Computer Simulation 37 

Appendix II - Graphs of Simulated System 
Trajectories 44 

Bibliography 68 



LIST OF ILLUSTRATIONS 



Pag© 



Figure 

(l-l) Geometric Interpretation of the Various Types of 
Stability 

(2-1) Control Function Uj As an Ideal Relay 

(2-2) Control Function U^ as a Relay With Dead Zone 

(2-3) Plot of y^(t) vs. Time For = -a^ sgn 

Stable Eigenvalue 

(2-4) Plot of y^(t) vs. Time For - -a^ sgn 

Unstable Eigenvalue, yi(0) > | & i j 

(2-5) Plot of y*(t) vs. Time For u^ - -a^ sgn 

Unstable Eigenvalue, y^(0) < | a i j 

(2-6) Control Function u^ With the Characteristics of a 
Saturating Amplifier 

(2-7) Geometric Interpretation of Successive Linear 
Approximation Process 

(4-1) Phase Portrait For an Uncontrolled VanderPol Equation 

(4-2) Root Locus Plot for Equation (4-2), - 00 < k < 0° 

(4-3) Root Locus of the Uncontrolled VanderFoi Equation in 
Stable Limit Cycle 

(4-4) System Optimum Switching Line 

(4-5) Root Locus for Uncontrolled System, Initial Conditions 

x — .5, x s .5 

Root Locus for Controlled System, Initial Conditions 
x = .5, x - .5 



1 

9 

9 

10 

11 

11 

12 

14 

20 

21 

22 

33 

35 

35 



(4-6) 



lo Introduction 



The stability of dynamic systems whose equations of motion are such 
that their solutions are difficult or impossible to obtain in closed form 
may successfully be analyzed by means of Lyapunov's Second Methods 

We have a given dynamic system described by the set of differential 
equations as 



x * i,( x ); jl(°) s 0 ° (1-1) 

Then* referring to Fig, (l-l), the origins, according to LaSalle and 
Lefschetz® is? 

stable whenever for each R there is an r «fR such 
that if a path (a motion) g* initiates at a point x° of 
the spherical region S(r) then it remains in the spherical 
region S(R) ever afters that is, a path starting in S(r) 
never reaches the bouncfeiy sphere H(R) of S(R)s 

asymptotically stable whenever it is stable and ia addition 
every path g^ starting inside some S(R C ), R 0 ^? o, tends to the 
origin as time increases indefinitely! 

unstable whenever for some R and any r, no matter how small, 
there is always in the spherical region S(r) a point x 
such that the path g - * through x reaches the boundary sphere 
H(R) e 




Figo (l-l ) 5 



Geometric Interpretation of the Various Types 
of Stability, 
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Now suppose we have a given dynamic system described by the set of 



differential equations written in vector form as 

x = £(x,u) 



( 1 - 2 ) 



or 

x = Fx + Du + £ (x,u) (1-5) 

where F is a constant square matrix, D is the control distribution 

* 

matrix, ^ is the control function vector, and £ (x,u) contains terms of 
a second and higher power in x and u. The linearized form of the equa- 
tion is 

x = Ex + Du (1-4) 

and can be said to adequately describe the motion of the system for 
small values of x. 

According to Kalman and Bertram^, a Lyapunov function may be de- 
fined as some scalar function of the state variables V(x) with the fol- 
lowing properties: (a) V(x) ) 0, V(x)( 0 when x £ x 0 , where Xg is the 

o 

equilibrium state, and (b) V(x) = V(x) s 0 when x - x q C 

We wish to apply Lyapunov' s Second Method to the design of a con- 
troller so as to return the state vector (x) to a position of equilib- 
rium (Xg) in some ootimum manner. No loss of generality will result 
if we choose (xg) to be the origin of the state space, and henceforth 
the equilibrium position will be considered as such. 

Lyapunov stated that if a Lyapunov funotion V(x) exists (being 
positive definite) in some neighborhood of the origin, and if -V(x) 
is also positive definite in some neighborhood of the origin, the 
system is asymptotically stable^. 

In general the met hod employed will be to make V(x) as negative as 
possible in order to dri ve the system to the origin as quickly as 
possible. 
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The Second Method of Lyapunov leaves open to the individual inves- 
tigator the choice of a particular Lyapunov function to be used in the 
analysis of a given system* For a given system then* there are possibly 
an unlimited number of such functionso The quadratic form of the canoni- 
cal state variables as the particular Lyapunov function employed in 
the design of the controller was governed by the simplification obtained 
in the control logic* 

Analysis and design for a controller of a linear , stationary system 
using the quadratic form procedes in a straightf orward manner 0 In the 
case of a non-linear system however, the limitation that the linear ap- 
proximation is sufficiently accurate only within some region surrounding 
the equilibrium position poses a rather severe restriction on the extent 
of the state space for which a given control can be useful* 

In order to overcome this limitation to some degree the controller 
was designed in such a manner that the control logic was obtained by 
successive linear approximations to the actual system at points on the 
system trajectory# 

For best performance of the system the time involved in the calcula- 
tion of the control logic should be infinitesimal so that there would be 
no appreciable time delay in the application of the proper control func- 
tion* That is, the system state vector would not have moved from the 
point (x^), where the state variables were sampled, to a new position on 
the trajectory, (x£) where the control function was actually applied 0 

In actual practice a finite calculation time is of course necessary 
so that the state vector will always have moved from the point of sampling 
before the new value of control function can be applied# 

For plants with long time constants, such as might exist in chemical 
or metallurgical processes or ship stabilization systems, presently 
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available digital computers have computation speeds adequate to success*, 
fully provide on line, real time control. 

Application of this method to plants with short time constants will 
depend upon increased speeds of the computers and the development of more 
refined methods of calculation. 
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2. Theory. 

If we have given a Lyapunov function V(x) s x^Px we may assert 
that^ the equilibrium state x e » 0 of a linear dynamic system 

x * Fbc (2-1) 

is asymptotically stable if and only if given any symmetric positive- 
definite matrix Q, there exists a symmetric, positive definite matrix P 
which is the unique solution of the set of linear equations described by 

pTp + pp = — Q . (2-2) 



A simple procedure for calculating the control logic attributed to 
R.W. Bass^ procedes as follows: Choose Q^O arbitrarily. P>0 may then 

be calculated from (2-2). ||x||2p will then be a Lyapunov function for the 

system (2-1). We now choose u(t) so as to make V (x) as negative as 
possible. Let 

V(x) » x^Px = ||x|| 2 P, then 

V(x) a X^Px + x^Px • (2-3) 

We wish to consider the effect of control on V. Taking the transpose 
of (1-4) we obtain 

x T » x T F T + u^D 1 (2-4) 



and 



V(x) ~ x^F^Px + x^PFx + u^D^Px + x^PDu 



or 

V(x) = x T (F T P + PF)x 4 2u T D T Px . 

T 

Since P is symmetric, DP * PD. Simplifying this gives 

V(x) » -x^Qx + 2u_^D^Px 
or V(x) * “|| X ||^Q ♦ 2u^D^ftc . 

To make V(x) as negative as possible we want the term 2u D Px to be 
as negative as possible. Therefore the control logic is given as 



(2-5) 
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u^(t) s -a^ sgn (D^Px)^ (2-6) 

where the are constants denoting the maximum allowable magnitudes of 

each individual control effort u^. It is obvious from (2-6) that in 

T T 

order to make the term Zu^D^Px as negative as possible the individual 
controller u^(t) must always exert a maximu m effort , and in such a 
direction that the term 



-a* 



sgn 



(D T Px) 



±'i 



is, in fact, negative. Thus the control logic becomes merely logic for 
detecting the appropriate switching points of the "bang-bang" controller 
Ui(t). 

In the case where there is only one controller the vector _u may be 
of the form such that uj * (0,0,0 . • . u). The system equation is 

then of the form 



(2-7) 



~ 0 1 0 ... 0 
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where the b^ are the coefficients of the characteristic equation of the 
unforced system, taken in reverse order. 

If we make the equivalence transformation from the physical state 
space to a canonical state space denoted by the matrix equation 

Z = Ox ( 2 - 8 ) 

then by differentiating we get 

Z 3 G * (2-9) 

and equation (l-4) becomes 

£ = G^FGz + G _1 Du . 
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We choose G such that 



G _1 FG = J\. 



( 2 - 10 ) 



where -A. is a diagonal matrix with the same eigenvalues as 

For the system with one controller such as described by (2-7) we 
may further choose the particular transformation matrix G such that 

G* 1 Du = E 



where s (a, a, a, . . .a). This has the effect of weighting the 

value of all the state variables equally in determining the sign of u(t). 

If the P of (2-2) is chosen to be the identity matrix I and the «A= 
matrix substituted for F, equation (2-2) becomes 

X T ! + 1 A. « -Q. (2-11) 



If all the real parts of the eigenvalues of F are negative, so then are 
all those of _A_ , and Q is positive definite. Also 

$(x) = £ T (_A_ T + -A, )%_ -fr 2E T £ . (2-12) 

The first term on the right side of (2-12) is again negative and the 
control logic is given by the second term. For a single controller, 
we have 

n 

u(t ) - -a sgn ( £ y±) • (2-13) 

In other words the switching function is such that the control must 
change sign whenever the sum of the n canonical state variables changes 
sign. 



Consider the given form 

“ -A.Z + g ” 1d ^ • 



#It is always possible to perform the equivalence transformation (2-10) 
provided F has no repeated eigenvalues. A system described by an nth order 
differential equation having repeated roots may be transformed into a par- 
tially uncoupled set of 1st order differential equations. It will be 
assumed that the systems considered in this paper have no multiple roots. 
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If we take the LaPlace transform we get 

sY(s) - Z (0) = _A_Y(s) + G“^-D^(s) 

or 

Y(s)(sl - J^) = ^(0) * GT ^D_u( s ) . (2-14) 

Since from (2-6) the control variable should assume only values of + a 

or - ^ , we stipulate that it is a constant during a given control 

period (between switchings ). We then have the LaPlace identity 

( G“ lD_u) ( s ) — _Z_ • 

s 

Solving (2-14) for a given component y^ we get 

Yi(s)(s -Xj.) - yi(o) + — — 



or 



Yi(s) = 



Yi(0) 



(s -\ ± ) 

Then the time domain solution of (2-15) is 

y*(t) - yi(0)e*^ it - 



s(s -X.) 



(2-15) 



*i 



(1 -.*1* 



) 



or 

Yi^) = -fyi(°) + --i- l e^ 1 - — . (2-16) 

*i J X ^ 

An interesting result becomes apparent from this solution as t in- 
creases without bound. For the case where X^ is negative (a consequence 

of Q being positive definite) we see that 

, . a. 

y i (t) - - } t oo 

x i 

This result indicates that although this form of controller is useful in 
making 'O'(x) more negative in an asymptotically stable system, it will not 
result in V(x) -*> 0 as t -► OO if u is a function of the type seen in 



Fig. (2-1). 
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a x' 


Control Function 


7- l U J 1, 


As An Ideal Relay. 



In fact it will result in a condition of chatter-motion around the origin 
as the controller becomes dominant at points in the state space where the 
yi<||e||s where ||e|| assumes some particular small value. 

This objectionable feature may be overcome in the case of asymptoti- 
cally stable systems by choosing the control function of the form shown 
in Fig. (2-2) below so that the controller is effectively disengaged 
within some region arbitrarily close to the origin, allowing the system 
to reach the equilibrium condition in an unforced manner. 




Fig. (2-2) Control Function As A Relay With Dead Zone. 

If the uncontrolled system were such that one (or more) roots con- 
tained positive real parts, (in this case -Q can no longer be positive 
definite) then for some particular points in state space, it is still 
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possible that 



v(x) = -||x 2 ||q . 

However V(x)>0 and remains so* 

Consider a plant governed by the equations 
dx/dt = Fx + Du(t) 

where the control variables are subject to the constraints, 

|ui(t)| < < co (i s 1, . . m) . (2=17) 

This is essentially the relay or saturating servo problem* 

The problem is to return every initial state to the origin* 

It is well known that if F has eigenvalues with posi- 
tive real Darts then there are some states which cannot be 
returned to the origin by any control subject to the con- 
straints (2-17) 4 . 

If a controller of the form described is applied to the unstable 
system, it is apparent from (2-5) that control may be maintained so long 
as the inequality (2-18) below exists. 2 That is 

-||x|| 2 Q + 2u t D t Px < 0. (2-18) 

From the solutions (2-16) the behavior of each y^(t) for 0<t<k, 
where k denotes the switching period is as in Fig* (2-3). 




Fig. (2-3) Plot of y^(t) vs. Time For u^ = -a^ sgn (jf^i* 

Stable Eigenvalue. 
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i 



a; 



\ 



Fig. (2-4) Plot of y^(t) vs. Time For * -a^ sgn (E^)^ 

Unstable Eigenvalue, y-; (0) > | a i |. 

I Xi I 




m 

Fig. (2-5) Plot of y^(t) vs. Time For = -a i sgn (E.^)^ 
Unstable Eigenvalue, yj(0) <| a i 1 . 

I X i I 

In the event that the system distribution matrix D were such that 
the controls were uncoupled, i.e. for each y^ there were a corresponding 
u^, then it would be possible to design a controller whose individual 
elements changed sign at the zero crossings of the individual state 
variables. 
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From Figo (2-5) it can be seen that the element would switch at 
t s To It is possible to calculate * ^(y-j^O)) from (2-16) for 

the case where y^(0)>0 # <0 . 

yi (t) s J y i (0) - !^L=le^ iT * — ii- , let y A (t) s 0 

— i Jr i 

then 



giving 








Xi 

■ ■ iT — ■ ■ — — . 



and 








(2-19) 



From (2-19) it is clear that for y^(0) * 0* T s — _ — . ln(l) - 0 

indicating that this system also will result in chatter-motion around 
the origin. 

Selection of a control function of the form of Fig, (2-2) would 
be unsatisfactory in this case because the uncontrolled system is un- 
stable at the origin, A better form of function would be that of the 
form of a saturating amplifier with a high gain 0 (See Fig, (2-6)). 

A u i ' 

- (ins). 

A> 

! 

Fig, (2-6) Control Function With the Characteristics of a 
Saturating Amplifier. 
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Employing LaSalle and Lefschetz's definitions^ of stability des- 
cribed in Section 1», a system whose eigenvalues all have negative real 
parts is asymptotically stable with no controls merely stable with a 
"bang-bang" controller (Figo (2-l)), and again asymptotically stable if 
the controller contains a dead zone as in Fig# (2-2) o Given an arbit- 
rary H(R), a system with eigenvalues lying in the right half plane is 
defined as unstable, and the same system with a controller of the form 
described may be defined as stable if X° is restricted to lie within 
some S(r) whose radius r is of some necessarily small value# The 
restrictions on the allowable values of r are necessarily related to 
•the constraints on the magnitudes of the control functions (2-17)» 

Given a non-linear system described by the equation x * F(x)x, 
the linear approximation of the non-linear system referred to by Kalman 
and Bertram may be obtained by a Taylor series expansion of the coef- 
ficient matrix F(x) at the origin# The linearized, constant, matrix F 

^ P „ / 

has the elements — i = « (The matrix F is the Jacobian of F(x) 

o Xf 

with respect to x at x s x(0) * 0) « 

Suppose a system described by the equation 
x = F(x)x + Du 

.is to be controlled in an environment where the disturbances are of 
such a magnitude as to carry the state variable beyond the region of 
validity of the linear approximation given by (l-4)# Referring to 
Fig# (2-7), this condition corresponds to an x° lying outside of S(r) 0 
By expanding F(x) about the point x° and obtaining the Jacobian of 
F(x) with respect to x at x s x° we may obtain a linear approxima- 
tion at x° 



c ^ 

x s F x ♦ Du 



( 2 - 20 ) 
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Figo (2-7) Geometric Interpretation, of Successive Linear 
Approximation Process « 

Control logic based on this approximation will be valid while the 
trajectory remains within the spherical region S(p)<» If S(p n ) is the 
region of validity of the linear approximation due to the n^* 1 sampling p 
then by an iterative process such that the trajectory from X a to X n+ “ 
always remains inside S(p a ) we may obtain an optimum control policy 
based on the particular Lyapunov function employed for the complete 
trajectory for x° to x e * 0 o 
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3 e Application of Theory to Design of a Single Controller* 

It was shown in Section 2* that the control logic resulting from 

the choice of the quadratic form of the uncoupled state variables as 

the system Lyapunov function is of the form 

u(t ) « -a sgn | i ^p 1 y* 

when the system has only one controller® Thus, the control problem is 

reduced to determining the sign of the sum of canonical state variables 

and insuring that the control, u(t), assumes the opposite sign.. The 

steps involved in this process consist of the followings 

I® Measuring the physical state variables (x)« 

2« Calculating the canonical state variables {%) from the 
measured (x). 

3. Forming ^y^ • 

4. Applying the proper control u(t) ^ -a sgn^yj.® 

As was shown in Section 2., the variables (^) may be calculated by 
the equation 

£ s • (3-1) 

-1 

The problem then reduces to evaluating G , 

It may be shown that a matrix G such that 

G" 1 FG = _A_ (3-2) 

v<here_/\-i s a diagonal matrix, may be formed by calculating the eigen- 
vectors of the matrix F and constructing an array in which each eigen- 
vector is a particular column of the array# Once a G has been obtained 
it is merely necessary to calculate the inverse, G~~, in order to be able 
to calculate the state variables (jr)# 

The elements of the diagonal matrix j\= are the eigenvalues of F# 

If the system is oscillatory the elements of _A. will necessarily be 
complex. Since the elements of the coefficient matrix F of the original 
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differential equation are all real, consideration of (3=2) reveals that 
for an oscillatory system there must be complex elements in G and G“^ s 

When G is complex, care must be taken in evaluating G“^ t If we 
let G a (feG + j$G, and attempt to form G* 1 - (&G)- 1 ♦ j(<3. G)* 1 

difficulty arises immediately 0 Since complex roots arise in conjugate 
pairs, there will be two eigenvectors (columns) of (t&G) due to the 
identical real parts of the conjugate pairso These vectors will be 
identical, and therefore linearly related so that (CfeG) will be singulars 
Similar reasoning shows that (c$ G) will also be singular, so that it is 
impossible to obtain (&G)“'^ and (<D G)~^ a 

Let us form a new matrix H in the f ollovdng manner e If the order 
of G is n, the order of H will be 2n so that for every complex element 
of G there will be four elements of H arranged in a square pattern© The 
two elements on the principal diagonal of the square have a value equal 
to the real part of the corresponding complex element of Go The ab- 
solute values of the remaining two elements are equal to the imaginary 
part of the corresponding complex element of G, the element in the lower 
left hand corner, taking the positive sign, while that in the upper right 
hand corner, the negative sign e Where the corresponding element in G 
is real, zeros appear in the locations in H corresponding to the missing 
imaginary parts 0 

If now we form then we may obtain the elements of from 

the proper locations of 

However since we are interested in obtaining the (y) corresponding 
to the measured (x), the simplest procedure is to form the 2 x a sup- 
plementary vector (x) which has alternately the real and imaginary 
parts of the physical state variableso Since these state variables may 
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be measured and have physical significance, they may only be real, so 

♦ 

that (x) consists of the values of the state variables alternating with 
zeros. Matrix multiplication of H” 4 on the right by ( x ) results in a 
2 x n column vector whose elements are the real and imaginary parts of 
the canonical space variables (%) , 

EXAMPLE , The product of a complex 2x2 square matrix and a 
column vector with real elements 



n 

yz 



ii 



+ o 1 
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l 12 4 ^12 



l 21 + J I 21 **22 4 



22 



x 2 



Then 



yi 

yz 

Rearranging gives 



R ll x l 4 ^ll^l 4 R 12 x 2 4 J I 12 X 2 

R 21 x 1 + jI 21 X l 4 R 22 X 2 4 jI 22 X 2 



yi - 


( R n x i + 


r 12 x 2 ) 


4 J( I ll 3C l 4 I 12 x 2^ 


( 3=3 ) 


y2 3 


(% x l 4 


^ZZ X Z^ 


4 *j^ I 21 X l 4 J 22 x 2 ^ 


(3»4) 



Suppose instead we form a 4 x 4 matrix of the form H"*, Then 



flfcyj, 




R n 


*'j I ll 


R 12 


-^12 




x i 


JLy x 




4 ^n 


R 11 




r 12 


* 


0 


&y 2 




R 21 


“J I 21 


R 22 


-jl 2 2 


x 2 


<*y 2 




3^21 


% 


*^ I 22 


R 22 




L °. 



and 

& Vl - R n x 1 4 R 12 x 2 

Ayi s ^ll^l. 4 ^1Z X 2 

ay 2 s % x l 4 R 22 x 2 

* A y 2 ■ V 21*1 4 j I 22 x 2 • 

Rearranging gives 

yi * ( R ll x l 4 R 12 x 2^ 4 j^ll x l 4 I 12 x 2^ 
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V2 * ( R 21 x l + R 22 x 2^ + ^ I 21 x l + I 22 3C 2^ 

which is identical to (3-3) and (3-4). 

It is evident from the above discussion that in the case of an oscil- 
latory system the transformation 

£ * (3-5) 

may map a given physioal state variable x^ from a point on the real 
axis to some point in the complex plane. Thus for a given coordinate 
axis in the physical n-space there is a complex plane in the canonical 
spaoe. 

Since, however, there is a one for one correspondence between points 
in the two spaces, a control which succeeds in reducing all state vari- 
ables to zero in the canonical space will have also reduced the real 
state variables to zero. 

The computation method employed to determine the transformation 
matrix G was based on a computer program MA.TSUB, programmed by 
Louis W. Ehrlich of the University of Texas, which calculates eigen- 
values and eigenvectors of a matrix up to order 50. The general method 
employed is to make an initial guess for each eigenvector of the matrix 
P and to converge on the correct value by an iterative scheme due to 
E.E. Osborne . The original program employed a guess of 1.0 for all com- 
ponents of the eigenvectors. Since it was assumed that the system 
matrix F was not changing at more than a moderate rate, a modification 
was made to the program so that the most recently calculated value for 
each eigenvector was selected as the initial guess for a new calculation 
thereby speeding the computation. The output of the MATSUB program was 
also modified so that the eigenvectors were arranged into two n x n 
matrices corresponding to the real and imaginary parts of the transfor- 



mation matrix G 
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From the two n x n matrices the 2n x 2n matrix H was formed and the 



inverse taken. The canonical state variables were then determined and 
the proper control calculated and applied. 

The path of the system trajectory in the physical space (x) was 
then calculated for a selected time interval between sampling instants 
by the Runge-Kutta-Gill method. Instantaneous values of the elements 
of the coefficient matrix F = F(x) were calculated for the new sam- 
pling point and the procedure repeated until the trajectory in the 
canonical space was within a pre-set epsilon region of the origin or 
a given time had elapsed. 

Initially the ideal case was simulated by holding the solution of 
the system trajectory at the sampling point until the new calculation 
of the control function was made and the control applied. Simulation of 
the non-ideal case, a finite calculation time, was affected by applying 
the control calculated at the (n-l)^ sampling, at the time of the n^b 
sampling. This is equivalent to a controller continually applying a 
new control as soon as it can perform the sampling and control cal- 
culations . 

Simplified flow charts of both the ideal and non-ideal case simula- 
tions are presented in Appendix I. The computer program in FORTRAN 
language (for the ideal case) is also included in Appendix I. 
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4. Simulation of Control of a Non-Linear Second Order System. 

A. Background. 

The non-linear VanderPol equation 

X - (1 -x 2 ) x + x * 0 (4-1) 

was selected as the system on which to evaluate the controller. Al- 
though it is only a second order equation it was felt that since it 
is quite non-linear it would be a fair test of the controller and the 
two dimension phase space has the advantage of making graphical analysis 
more simple. 

The equation exhibits a stable limit cycle of the general shape 
shown in Fig. (4-1). • 




Fig. (4-1) Phase Portrait For an Uncontrolled VanderPol Equation. 
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; 




Referring to (4-l), when x ■ 0, the system is most unstable. For 
this point the roots are +.5 ± j.866. It was determined from trajec- 
tories of the uncontrolled system obtained by digital computer analysis 
that when the system was in the limit cycle the maximum excursion of 
displacement was ~ 2.1925. This value was substantiated by an analog 

computer simulation. From (4-l) the eigenvalues, or roots, of the system 
vrhen the displacement reaches its maximum value in the limit cycle were 
calculated to be -3.5235 and -.2835 so that the locus of roots of equa- 
tion (4-1 ) during a complete circuit of its stable limit cycle is shown 
in Fig. (4-3). . • 




Fig. (4-3) Root Locus of the Uncontrolled VanderPol Equation 
in Stable Limit Cycle. 

Examination of Fig. (4-l) or equation (4-l) reveals that the system be- 
comes unstable whenever |x| < 1.0 and becomes stable again whenever 
|x|> 1.0. Thus the equation presents the controller with the problem 
of controlling a system which contains at various times, stable real 
roots, stable complex roots, and unstable complex roots. 

It is informative to observe that a stability analysis of the 
VanderPol equation by Lyapunov's Second Method yields precisely the 
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same results as the more familiar methods of analysis such as the root 
locus. 



Rewrite (4-1) in vector form. Then 



or 



and 



Choose 



then 



i • [-? <^>] z. 



*1 1 *2 



X 2 = -X]^ + (l-x^)x. 



2 * 



V(x) = j|x|| • xj + x§ 



V(x) ■ 2x^xi + 2x 2 xg 
and from (4-3) and (4-4), 

'O'(x) = 2 x^x 2 + 2x 2 (-x^ + (1-Xj)x 2 ) 



(4-2) 

(4-3) 

(4-4) 



(4-5) 



or 

V(x) = 2x| ( 1-x^ ) (4-6) 

so that V(x) is positive definite but -V(x) is not for |x^j<l. There- 
fore the region near the origin (-l<x<l) is unstable. This is pre- 
cisely the result obtained by observing that the root locus passes into 
the right half plane whenever |xjJ <1. 

B. Procedure. 

It was decided to initially investigate the action of the control 
with no time delay between calculation of the control and its application 
in order to verify that the method of control and the simulation pro- 
cedure (described in Section 3.) were feasible. First a family of tra- 
jectories of the system with no control effort was obtained from 
various initial conditions chosen inside and outside the limit cycle 
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(Graphs A-l through A-8). Trajectories were then obtained using the same 
initial conditions with the control system activated (Graphs B-l through 
B-8). The sampling rate was set at .3 seconds and u was set at 1.0* An 
arbitrary epsilon region surrounding the origin was established so that 
the trajectory was considered to have reached the origin when the sum 
of the canonical state variables squared was less than .001. Whenever 
the trajectory entered the epsilon region of the origin the solution 
was terminated. 

Effects on the trajectories due to varying sampling rates (Graphs 
C-l through C-5) were then investigated, and finally two trajectories were 
obtained for the system simulating the effect of finite calculation time 
(Graphs D-l and D-2). These trajectories all were started from an 
initial condition close to the uncontrolled limit cycle trajectory. 

The simulation program included a sub-rort ine for the graphical pre- 
sentation of the system trajectories. The graphs noted above are pre- 
sented in Appendix II. 
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UNCONTROLLED CASE: SAMPLING RATS s 0.3 SECOND. 
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II-CU CONTROLLED CASE: SAMPLING RATE a .50 SECONDS, u = 1*0 
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TIME-DELAY CASE: INITIAL CONDITIONS x ' = 1.5, X = 1.5, u a 1,0 
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D. Discussion of Digital Computer Simulation , 

Comparison of the trajectories of the controlled and uncontrolled 
systems starting from the same initial conditions (Graphs A-l through 
A-8 and B-l through B-8) clearly indicates that the control system with 
a control effort = 1.0 and a sampling rate of .30 seconds^' is capable 
of driving the system into the unstable region surrounding the origin and 
maintaining it within some region of lesser extent than that enclosed by 
the system’s uncontrolled limit cycle. Referring to Pig. (l-l), if the 
radius of S(R) were established for this system as R * .5, then the con- 
trolled system might be defined as stable while the uncontrolled system 
would be defined as unstable. Graph B-l, initial conditions x ■ .2, 
x » .2, illustrates the system trajectory in the region close to the 
origin. The chatter-motion mode of operation is readily apparent. 

The theoretical analysis of Section 2., predicting that by employ- 
ment of this type of controller, an unstable system could be made stable, 
but not asymptotically stable, and that a chatter -mot ion mode of 
operation would result around the origin, is substantiated. 

Examination of Graphs C-l through C-5 reveals that an increase in the 
control system sampling rate causes the trajectory, having once arrived in 
the neighborhood of the. origin, to remain thereafter within a smaller re- 
gion of the origin than the same system operating with a low sampling rate. 

Curves B-3 and C-3, both trajectories for the identical initial condi- 
tions and system parameters, demonstrate remarkably different trajectories, 
each of which eventually arrives in a neighborhood close to the origin. 

#The time interval between successive calculations of the control 
function. The Runge-Kutta integration of the trajectory employed an 
interval of 0.03 seconds. 
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Examination of the calculated data revealed that in the case of 



Graph B-3 the initial values of the canonical state variables were 

= .75 - j.1561 

y2 = .75 + j.1561 

while in the case of Graph C-3 the values were 

y-L ■ .75 -j.1561 

Y2 = .75 —j.1561 « 

Since the control function u * -1.0 sgn ( E&yi + E<!yib 

it is evident that u assumed a different sign in each case. That this 
was so can be determined by examination of the diverging paths of the 
trajectories in the two cases. 

It seemed evident that the more direct trajectory should be the 
correct one, however the possibility was raised that there might not 
be a unique solution for the control function. Six additional runs 
each exactly coincided with the more direct trajectory B-3, indicating 
that oerhaps an error in calculation led to the selection of the op- 
posite value of control function rather than there being two possible 
solutions for the control function. 

The computation routine had been programmed so that the inverse 
transformation matrix was printed at each sampling point. It was dis- 
covered that the transformation matrices in the two cases were 
not identical. 



INITIAL INVERSE TRANSFORMATION MATRIX FOR GRAPH B-3 



.35820E-10 
-.64051E-00 
- .35820E-10 
.640516-00 



.64051E-00 
.35820E-10 
-.64051E-00 
- .35820E-10 



. 50000E-00 
-.40032E-00 
.50000E-00 
.40032E-00 



.40032E-00 
. 50000E-00 
-.40032E-00 
.50000E-00 
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INITIAL INVERSE TRANSFORMATION MATRIX FOR GRAPH C-3 



.37180E-08 .64051E-00 .50000E-00 .40032E-00 

-.64051E-00 .37180E-08 -.40032E-00 .50000E-00 

© 50000E-00 .40032E-00 .65670E-09 .64051S-00 

- .40032E-00 .50000E-00 -.64051E-00 .65484E-09 

In order to determine if two solutions were possible, the two in- 
verse matrices were checked against the original matrix to see if both 
were valid o 

The original matrix was calculated by hand in the same manner that 
the computer was programmed to do 0 

The procedure is illustrated below. 



The initial F matrix was 



fo 1 1 

L-l -1.25J 



and the initial eigenvalues were 



-.625 * j. 7806. Then for the first eigenvalue 

s (-.625 * o„7806) 



p [:;] 






or 







x 2 




(-.625 


* j.7806)x^ 






X 1 


* 


(-.625 


- o » 7806 ) x 2 


so that 
















x 2 




1.0 








X 1 


= 


-.625 


- o » 78 0 6 o 


In a similar 


manner, for 


the second 


eigenvalue 






x 2 


S3 


1.0 








X 1 


S 


-.625 


^ 0 • 7 606 , 


Then the 2n : 


k 2n matrix becomes 








-.625 




-.7806 


-.625 .7806 






.7806 




-.625 


-.7806 -.625 






1 




0 


1 0 






0 




1 


0 1 



When this matrix is multiplied by the inverses, B-3 and C-3, only the 
inverse B-3 results in the identity matrix, demonstrating that the 
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calculation resulting in C-3 is not an independent correct solution* 

Comparison of the false trajectory with the correct one indicated 
the possibility of there being at least two other instances of erroneous 
calculations besides the one at the origin* These are indicated by the 
abrupt changes in direction (outward) of the trajectory at the points 
x - .65, y s 0,6, and x £ '.5, y s *45 . Calculations similar 
to the one carried out above verified the fact that erroneous determina- 
tions of the correct control function had been made* 

The similarity of the trajectories of curves C-l, C-2 and D«X in 
the fourth quadrant raises the interesting comparison between the control 
obtained by Lyapunov's method and the method of optimum control employing 
the calculation of the optimum switching lines in negative time as dis- 
cussed by I. Flugge-Lotz and H.A. Titus^. 

We may speculate that there is some optimum switching line for 
this system, emanating from the origin in a manner similar to that 
in Fig* (4-4). * 




Fig. (4-4) System Optimum Switching Line. 
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The optimum trajectory for this system would then proceed from the 
initial condition to its first intersection with the optimum switching 
line, change the sign of its controller, and proceed along the optimum 
switching line to the origin. 

It appears evident that the method of switching logic discussed in 
this paper selects a switching time too soon in the cases of curves C~1 
and C-2, In the case of curve D-2 the system appeared to be in a state 
of indecision. The switchings at points A and C undoubtedly would have 
brought the trajectory initially closer to the origin than the one 
finally selected at E. The switchings at B and D served only to hinder 
the system response. Investigation of the computer data again indicated 
erroneous calculations made for switchings at points B, C and D, although 
the resulting decision at C proved to be a correct one. 

Graphs D-l and D-2, showing the simulation of the system operating 
under non-ideal conditions, again substantiate the fact that higher 
sampling rates generally produce more satisfactory operation of the 
system. It is interesting to note that the controlled system operating 
with a time delay of .3 seconds eventually settles into a stable limit 
cycle of approximately one half the size of the uncontrolled system. 

The effect of the controller on the excursions of the system roots 
as the state point travels along a trajectory for a particular initial 
condition is demonstrated in Figs. (4-5) and (4-6). The numbers indicate 
time in seconds to reach the position indicated. 

The root locus for the uncontrolled system starting from initial 
conditions x = .5, x = .5 is illustrated in Fig. (4-5). The roots 

attain their most stable positions after the trajectory has settled into 
its stable limit cycle. 
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Fig, (4-5) Root Locus for Uncontrolled System, Initial 
Conditions x = ,5, x ~ ,5, 

As the controlled trajectory approaches the region of chatter- 
motion surrounding the origin, the root locations depart very little 
from their most unstable position (+.5 ± j«866) because maximum 

instability occurs when x ~ 0, The most stable root positions on 

the controlled system root locus occur early in the solution while 
the trajectory is relatively far from the origin. 




Fig, (4-6) Root Locus for Controlled System, Initial 
Conditions x = .5, x s ,5. 
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5. Conclusions* 



The operation of the controller in regulating a simulated non* 
linear system described by the VanderPol equation supports the theo- 
retical predictions developed in the design* The solution obtained by 
this method is seen to be less than optimum because the logic does not 
cause the control to switch at the optimum svdtching curve* It is 
apparent from the figures that the sampling rate for determining the 
eigenvalues, and hence changes in control action, is critical for this 
to be an effective control design procedure* With the fine sampling 
the system disturbances were rapidly and effectively controlled* 

There exist today very few direct techniques for the control of 
non-linear systems* The procedures studied here provide such a tech- 
nique* The method has the advantage of being applicable to non* 
linear time varying systems of any order* It does however require a 
digital computer in the control link* 
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APPENDIX I 

FLOW CHARTS AND FORTRAN LANGUAGE PROGRAM 
FOR DIGITAL COMPUTER SIMULATION 
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Fig. 1-1 FLOW CHART FOR IDEAL CASE SIMULATION. 
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Fig. 1-2 FLOW CHART FOR FINITE CALCULATION TIME SIMULATION 
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.23 

fHi 1 rnfUc ?ii?m ffimiu mwc^ 

SjofF^S^tF^S^P^^I^LSi^oS'?? C ° EFF MATRIX 

1204 F 0 R rAr!E 22 .'Ii C 3 E ! i , 2 N 5 ) RSe IRA «SFORHATi ON MATRIX//) 

1205 FORMAT (1 HO) 

7002 STOP 

1052 FpRKAi-?^’,Mo?"’ J1 ’ J = , ’ Km ’ I = , ’ MMI 
PRINT 1200 
L I N E = 4 
TO = 0.0 
NUMPTS =0 
LAG = 0 

DO 200 J = 1 . MM 
DO 200 I = 1 ,&M 
GR ( I , J ) = 1 . 0 • •; 

G I { 1 v J ) = 0 .Q; 

A I l I » J ) = 0.0 

T = TO , , 

!l ' EC0TRS 0F AR AN ° F0RHS 

/100 ° nn L kn 9 S V B -S M » ALRS »ALIS) 

D 0 = 30 J = 1 ,„« 

K = M -1 

. XEXP (K, 1 ) = XR II) 

3 o L xFxp N ?M'?i l > ir i N! o o ORMS D0UBLE SI2E transformation MATRIX CALLED DUB 

DO 31 1 = 1 , MM * 

M = 2*1 

K = M-l ✓ 

DO 31 J = 1 , MM 
L= 2 * J -1 
N= 2 * J 

DUB { ,t , L ) = GR ( I,J) 

DUB (M,N) = GR(I ,J) 

,, -DU 6 (K,N) =-G I ( I , J 5 

3 -' 2 G,U ' J ' 

c of E 5 ub , = ou dubiw I!inv ’ KEY1 

2 PRINT | , ' 2, ’ KtY 
STOP 

? ^RSr-J.U- 9H GR MA TRIX SINGULAR) • 

^ I UJN I t N U c 

F0R gg 32 ^I=UNN STATE VARIABLES YEXP (CONTAINS BOTH REAL AND IMAG. TERMS 
YE X P ( I ) = 0.0 



JOO 



40 



CO ^ 2 J — ] 7 N M 

32 YSQR { = ) C =: C YEXP(! } + OtJBINVII.Jl * X5 XP{ J , 1 ) 

00^1003 I = ] » \i M 

noo zbr-' uoo ’ ,40, ’ ,Eo ' ‘ 

PRINT' 1203 

LINP^’l^N^ { - D 9 BiX ' V(I ’ J),J = 1,NN5,I = 1 ’ NN5 

1402 pIuI? Taft 1 "* 

LINE = 4 
GO TO 1404 
14 03 PRIM 1205 

LINE = LINE 3 

c , , T£STS for closeness to origin 

] 005 ftlK* : J5CRI 1005, 1303, 1303 

C CALCULATION OF CONTROL FUNCTION 
DO 1006 I = 1 , NN 
1006 YRSUM = YRSUM + YEXP(I) 

UfiP = - S IGNF ' J , YRSUM ) 

IP if' tvALUES(2^ 1 5?fooT7oVV7o5 2 ’ 2 ’ >_,0 * E_6> 702 ' 702 ' 70 > 

700 VALUES ( 2 ? 1 1= -VAIUES(2,I) 

70! ?5fc¥f§il’ 2 *- - VALUESI2 ’ 2) 

TOEL = 7-«- DELTA 

1304 IF (ID EL - T - DT) 1300 » 1301 * 1301 

m isibsv.iev 302 '™-* 7 * 

GR l T0’l000 ^ * {U ° - XR{1) * XR(1) > i 

1303 CONTINUE ■' 

PRINT 21202 

21202 P0RMAT t {26hP-tRAJ£CT0RY INSIDE EPSILON///) 

* f\ l W I n J Z U i! 

h,2 ° 2 T ; YSCR R4DIUS VKT0R/// ’ 

G5APH'fNUNPTS,X!,Y!,8> 

c 1 5§f U Kfi^i5VW T f g« C ° N I R ? t !ii 0 , s 2f T f u5?f JECT0RV F0R PERIOD 

nos k\M ’“OS. '-06, 1406 

LINE = 4 

1406 PRINT 1 202 , T, XR ( 1 ) , XR { 2 ) 

LINE = LINE + 1 
I Ft 900 - NUXPT S 5 1304,750,750 
7b0 NUMPTS = NUN? T S + 1 ’ 

XII N'JMPTS) - X R { 1 ) 

Y! f NUMPTS) = XR { 2 ) 

CO f 0 1304 
END 

SUBROUTINE RKUT7A ( M, T , XR ,DT , UOP) 

CCNSopl 1 ,£?*ggi&f?iAt6ll' xt>OT « 2c> » xcl20,t Cl4 ’ 

C { 1 ) — 0 o 0 

c ; 2 > =0.5 

C t 3 ) = 0.5 
C ( 4 1 =1.0 
DO 1050 I =1,4 
TC = T + C ( I ) * OT 
DO 1051 J = 1 ? M 

*051 XC(J) = XR ( J ) •}• C ( I ) * AK { I — 1 ,j) 



CALL DERI V ( TC , XC , XDGT 



UOP 



) 



41 



+ 2.0* AK(2,J) + 2 . C* AK! 3 » J 5 +AK(4,J))/6. 

IT, UOP. M ) 

r AR! 20,20) ,AI (20,20 ),GRl 20,20 ),GI< 20, 



+ A R ( M , 2 ) * XR ( 2 ) + UOP 
ALIS ) 



DO 10 50 J =1,M 
1050 A !< ( I ■> J ) = DT * XDOT ( J ) 

DO 1052 J = 1 , M 
1052 XR(J) =XR ( J ) +(AK(1,J) 

END 

SUBROUTINE DERI V (T,XR, XD 
DIMENSION XDOT (20), XRI2C) 

10), VALUES! 2, 20) 

COMMON AR ,A1,GR,GI, VALUES 
XDOT! 1 ) = XR ( 2 ) 

XDOT! 2) = AR(H, 1 ) * XR(l) 

END 

SUBROUTINE GSUB (M , ALRS , ALii l „„ 

OD I MENS I ON AR( 20,20) .All 20,20 ) ,BR( 20, 20) ,BI (20,20) ,CR (20. 20) , 
1CI ! 20,20 ) ,XR( 20 ) ,XI (20 5 , YR ( 20) » YI ( 20 ) , VALUES ( 2 , 20 ) , GR ( 20 , 20 ) , 
2G I! 20 » 20 ) 

COMMON CR , Cl ,GR,GI .VALUES 
N = M 
LOOP = 1 

922 DO 519 I = 1 , N 
DO 519 J = 1 , N 
BRC I, J) = CR! I, J ) 

AR! I - J ) = CR! I , J ) 

= Cl ! I , J) 

519 AI ! I , J ) = C I ! I , J ) 

ASSIGN 917 TO IA 
ID 



811 TO 



TO ID 
530 TO IA 
40 TQ IB 
523 TO IC 



C 



ASSIGN 
MM=K 

GO TO 535 
GO TO ID 
ASSIGN 914 
ASSIGN 
ASSIGN 
ASSIGN 
l SL =- 1 

GO TO 92 ' ►. 

ISL = 0 
ALR = ALRS 
ALI = ALIS 
I T= 1 

__ DO 504 I = 1 , N 
NEXT TWO INSTRUCTIONS DESIGNATE INITIAL 
PRECEDING EIGENVECTOR SOLUTION 
XR! I )=GR! I jLOOP) 

XI { I)=GI ! I^OOP) 

DO 5 1 = 1 , M a 
AR! I , I ) = A R (I , I ) - ALR 
Alt I»I)=AI< I»I J-ALI 
I J = 1 



917 

811 



523 

9 



403 



504 

4 



GUESSES FOR EIGENVECTOR 



10 



1 1 



12 

13 

106 



1C9 



14 



B IG=0. 

DO 13 1 = 1, N 
YR! I )=0. 

YI ! I )=0. 

DO 11 J = 1 , N 
YR! I ) =YR ! I ) + AR ! I 
YI ! I )=YI ( I )+AI ( I - • 

AM = YR ( I )*YR( I )+YI ( I )*YI (1) 

IF (AM-BIG) 13,13,12 
8 1 G = AM 
J J= I 

CONTINUE 
IF (BIG) 109,106, 

I C T- 1 000 
JJ = N 
ISL = 1 
GO TO 92 
RQNR=0 . 

RQN 1=0 . 

R0D=0 . 

DO 14 1=1, N 
RQNR=RQNR + XR ! I ) *YR! I ) +X I ! I ) » Y I ( I ) 
RQN I = RQN I + XR ( I ) *YI ( I 5 -XI ( I )*YP. C I ) 
RQD=RQD+XR ( I) *XR ( I) +XI ( I ) - v 
AMUR=RQNR/RQD 



,J)*XR(J)-AI( I,J)*XI ! J) 
J)*XR( J)+AR{I,J)*XI (J) 



109 



-XI ( I ) 



L2 



AMU i -RQNI /RGD 
A t : M = A M U R - A M U R A M U I * A M U I 
81 TS = 0„ 

08 15 I = 1 , N 
150TS=TS- 

1 j y T { 

'r.0 1 l 6' 1 = 1 r N 

■' l M I >=< YR( JJ) e-YR( I5+YZC J. • ..YI { n )/Bi G 



1 6 

i 1.1 
13 
.19 

20 



XR 



5 I = 1 , N 

■ S + { YR.C I ) - AMUR* XR ( I )+AFUl*XI ( I 
I ) -AMUR* X I { I ) - AMU I »X R ( I } ) * *2 

6 I = 1 N 



) )* * 2 + 



fRi i n/ciG 



310 

99 

100 

29 

535 



21 



22 



XI (I )={ YR{ JJ)*YI { I )-Yi 
XRC JJ } = 1 .0 
XiiJJhC.C 

(TS / RQD,- 10.E-4) 20,20,18 
1 1- ( I J - 20) 19,20,20 

GO fn’io 
I C T= I J 

M I T 2 = IJ + 2 0 
ALR=AMUR+ALR 
All = AMU I 4-AL I 
'. V -N 

OG 310 1=1 ,N 
ARC I, I ) = AR { I , I ) - Af'lJR 
AI( I, I)=AI ( I , I )-AMUI 
GO TO 29 
DO ICO 1=1, N 
ARC I, I ) = AR t I r I ) - ALR 
Aid, I ) =AI (1,1 ) - A L I 
i J=I J-H 
00 27 1=2, MM 
I M 1 = 1-1 
DO 27 J = 1 , I M 1 

FM= A R C I , J ) * AR ( I , J 5 + AI ( I, J ) »AI { I , J ) 
S.:-AR( J, J 5* ARC J, J 5 + AI t J, J )*AI ( J. J ) 
r C F M- S M ) 24,24,22 
<0 23 K = J,MM 
T 1 = AR { J , K ) •. 

T 2 = A I ( J , K ) 

ARC J,K)=AR{ I , K) 

A I C J , K ) = A I ( I , K ) 

AR ( I , K ) = T 1 ' 

23 A I C I , K) =T2, 

T 1 = X R ( J ) 

T 2 = X I ( J ) 

XR £ J ) = X R ( I ) 

XI ( Ji=XI U) 

XR { I ) =T 1 ; 

XI { I ) =12 t 
T 1 = F M 

f;-’=sm 
sm= ri 

IF (SM) 25,27,25 
IF ( F M ) 90 » 27 - °0 

= JARC i,H! ; i ■< J ) +AIC I , J 5 * AI ( J , J ) ) /SM 

nn = i£ R A Jl 1 J A?, Al 1 I 5 J5-ARC I , J )*AI ( J , J ) ) / SM 
DO ^ o K=J«MM 

AR£I,K)=AR(I , K ) - RR* AR { J , l< } + R I * A I C J,K ) 

Al IfKJ-^nijKl-RRtA 1 (J,K)-RI»AR(J,K) 

AR ( I , j } =0 . 

AT C I , J ) =0* 

XR C I )=XR( I ) -RR*/R ( J )+RI*XI(J) 

XI { I 5=XI ( I )-RR*XI { J)-RI»XR( J) 



24 

25 
90 



26 

27 

550 



75 0 
75 1 

23 



GO TO IA 
SMALL = 1 000 * 

DO 28 K= 1 , MM 
I K K = K 

T2=AR(K,K)*AR(K,K)+AI (X,K)»AI(K,K) 
Ir (115 roO , 752, 750 
IF ( T 1 -SMALL 5 75 1,23,28 
S.M ALL = T 1 
I 0'-i< 

CONTINUE 
GO TO IB 



ko 



752 12 = IKK 

IF USD 753,30,30 

30 157=1 

!l i -2 000 

00 974 1=1, MM 
XR!I ) =0.0 

974 XI { 1 )=0.C 

753 YR( iZ 5 = 1 .0 
YI ( 12 5=0.0 
J J = I 2 

R 1 G= 1 . 0 

IF ( I Z - 1 * M ) 33,32,33 

32 122=2 

GO 70 95 

33 122=12+: 

DO 31 I = I Z 2 , MM 
YR( I)=0. 

31 YIU)=0. 

I ZZ=MM-I Z+2 

1 r U2-1) 95,49,95 

40 122=1 

41 •: IG-0. 

95 DO 4 6 I = I 2 2 , M M 
11= MM- I + 1 
XK = I I + 1 
SR = 0 . 
si=0 u 

IF (I-I) 42,44,42 

42 DO 43 K=KK,MM 

SR = SR+AR ( I 1 ,K)*YR(K)-AI ( I I , K ) * Y I (K) 

43 SI =SH-AR { I I ,K)*YI ( K) + AI { I I , K) *YR(K) 

44 T 1 =AR J 1 1 , 1 1) * AR { 1 1 , III +A I ( I 1 , 1 I ) *A I ( 1 1 , 1 1 1 

YR( 1 1 } = { AR( l l , I I ) *{ XR{ I I 3-SR1 +AI ( I I, I I ) *(XI ( I D-SI ) ) /T1 
YI { 1 1 } = { AR(I I v I I )*( XI ( I I 3 -S I >-AI ( I I, I I )* (XR( II )-SR ) ) /T1 
AM=YR ( I I ) »YR { 1 1 ) +YI ( I I ) « Y I U I) 

IF (AM-BIG) 46,46,45 

45 JJ=I! 

B I G=AM 

46 CONTINUE " 

49 DO 47 1=1, MM 

X R ( I )=( YR( JJ>*YR( I ) +YI ( JJ5*YI ( I ) )/8IG 

47 XI U}=(YR(4J)+YI (I)-YI ( JJ)*YR( I ) J/BIG 

xrj j j ) = i . 0 ; ' 

X I ( J J ) =0 . 0 >, * 

92 00 601 I = 1 , N 
DO 601 J = 1 , N 

AR ( I , J ) = B R ( I , J ) 

601 AI(I,J)=BI(I,J) 

116 IF USD 755,50,60 
755 GO TO IC 

50 ALR=0. 

A L I = 0 . 

S U M = 0 . 0 

55 DO 52 1=1, N 
YR( I )=0. 

YI ( I )=0. 

CO 51 K=1,N 

YR ( I 5 = YR ( I ) + AR( I ,K) *XR(K 5-AI { I ,K)*XI ( K ) 

51 YI II) =YI { I ) +AR { I , K ) * X I ( K ) + A I ( I , K ) * XR ( K ) 

ALR=ALR+XR( I )*YR( I ) + X I ( I ) *Y I ( I ) 

ALI = ALI+XRm«YI ( I ) -X I m«YR( I ) 

52 SUM=SUM+XR( I)«XR( I)+XI (I )«Xl (I) 

ALR=ALR/ SUM 

AU =AL 1/ SUM 
AM=ALR«ALR+ALI*ALI 
B3 TS=0 , 

DO 53 1=1 , N 

T I = V R ( I)-ALR*XR( l)+ALI*XI{ I ) 

T2 = Y I ( I )-Al.R*XI ( I )-ALI*XR( I ) 

53 TS=TS+T1*T1+T2*T2 

93 IF (IS / SUM - 10.E-14) 60,60,301 
301 IF (IJ-MI72) 99,400,400 

400 IF (IT - 3) 402,990,402 
402 ALR = ALR + .1 



60 

63 



All = All 
I T = I T + 1 
CO 10 4 
I S L = 0 



VALUES ( 1 ,LOOP) = ALR 
V A L U E S ( 2 , L 0 0 P ) = ALI 
IF { J J-M) 6] , 65, 61 
Ol T1=XR(JJ) 

T 2 = X I { JJ ) 

XR [ JJ ) = X R ( N ) 

X I f J J ) =x I { N!) 

X R ( N ) = T 1 
XI (N)=T2 
CO 63 K=1,N 
1 1 =AR { 

T 2 = A I { JJ, |< i 
AR{ JJ ,K ) =AR( M»K ) 

A I ( J J , K ) =A I { N , K ) 

AR ; N » K } = T ] 

63 A I l.\',KJ«72 

00 62 K = 1 j N 

1 3 — A R { K , J J ) 

! 2 = A I ( K ? 'J J ) 

AR(K» JJ)=AR(K,N) 

A I K, JJ)=AI (K,N) 

AR { K , N 3 = T1 
A I I K , N } = T2 
N=N- I 

DO 66 I = 1 , M 
DO 66 J = 1 , N 

= I ? J i)-XR{ I )*AR (N+l vJJ+XKI }*AIIN+1 n 

DO 600 J) ~ XR(I,<tAI (N+1.J)-XI(I )*AR(N+llj) 

00 600 J = 1 * Xi 
BR { I j J ) = AR{ I , J ) 
dlf It J)=AI ( I , J) 

DO 702 1 = 1 M 
DO A02 J = 1 1 M 



62 

65 



66 



600 

700 



ARC I,J)=CRfI,J) 



701 

702 



911 

914 

915 

703 



916 

704 



aic i ? j)=ci||; j) 



. ARfcjt; I 3-ALR' 
a1{I»I)=AI||- ; 



'TO IA 



I A 



if n-j) 

AR I I , I ) = 

P. » * . , I j 

CONTINUE 

M * ; .v; 

ASSIGN 911 
GO TO 535 
ASSIGN 530 TO 
I SL=- 1 

00 703 1=1 ,M 
XR ( I )=0. 

XI f I >=0. 

ASSIGN 753 
ASSIGN 7G4 
GO TO 530 
ASSIGN 525 
GO TO 92 



701.702 



-ALI 



TO 

TO 



IB 

IC 



TO IC 



'04 DO 920 I = 1,M 

EIGENVECTORS 'xR^ANO^'xi TRANSFORMATION MATRICES GR AND GI FROM 

GRUjLOOP) = XR ( I } 

Gilt, LOOP } = X[(l) 

ASSIGN 40 TO IB 
LOOP = LOOP + 1 
IfjtN-l) 921,67,523 
ALR=AR(1,1) 

AL I =A I (1 , 1 3 
VALUES ( 1 , LOOP ) = ALR 
n VALUES(2,LC0P) = ALI 

GO TO 700 
921 GO TO 4000 



920 



525 

67 



433 



XMAX = -10. «* 2 0 
03 1 000 I = J,M 

rn (XMAX - VALUE S ( 1 , I ) ) 3,3,1000 
3 X AX - VALUES! I, 15 

1 < X — I 

1CC0 COM! I VUE 

IF (MX-J) 1,4000,1 
i VALUES (1, MX) = VALUES! 1 , J ) 

V A L U E S ! 1 , J 5 = XMAX 
SAVE = VALUES (2, MX) 

VALUE S (2, MX) = VALUES ( 2, J) 

VALUES ( 2, J 5 = SAVE 
03 3000 K=1,M 
SAVE = OR (K, MX 5 
GfUK,MX3 = GR(KtJ) 

CR { K , J ) = SAVE 
SAVE = G I ( K , M X ) 

6I(K,HX)= 6I(K,J) 

3000 GI(K,J) = SAVE 
4000 CONTINUE 
GO TO 992 

990 PRINT 991 

991 FORMAT ( I 5H NO CONVERGENCE) 

992 END 

SUBROUTINE GAUSS3 ( N , E?,A, X,KER> 
DIMENSION A ( 40 , 40 5 , X!40,40) 

DO 1 1 = 1, N* 

DO 1 J = 1 , N 

1 X(!»J)=0.Q 
DO 2 is = 1 , N 

2 x { :< , ;< ) = U 0 

10 00 34 L = 1 9 N 
KP = G 

2 = 0.0 

DO 1 2 K= L , N 

IF(Z-ABSF(A(K,L) )) 1 1, 12, 12 

11 Z =A8SF ! A ! Kj-L ) 5 

j/ p — !^ ' r ■ 

12 CONTINUE > : 

IF(L-KP) 1 3 » 2 0 } 2 0 

13 CO 14 J=L,N 

Z = A ( L , J 5 " 

A ( L , J ) =A ( KP , J ) 

14 A ! K P , J ) = Z 

00 15 J= 1 » N 
Z = X ( L j J ) 

X ! L , J 3 =X i KP » J ) 

15 X ! KP f J } =Z 

20 IFtASSF! A!L,L ) ) -EP ) 50 , 50 , 30 

30 IF(L-N) 51,34,34 

31 L P 1 = L -J- 1 

DO 36 K= LP 1 , N 
I F { A { K , L ) ) 32,36,32 

32 RATIQ=A(K,L)/A(L,L) 

DO 33 J = l. P 1 , N 

33 A(K, J)=A(K, j)-RATIO*A(L» J) 

DO 35 J = I , N 

3 5 X(K, J)=XvK»J)-RATIO»X(L, J) 

36 CONTINUE 

34 CONTINUE 

40 Du 43 1 = 1 ,N 

1 I = N +1-1 

DO 43 J = I , N 
S = 0.0 

IF! 1 1 - M } 4 1 , 43,43 

41 I I P 1 = 1 1 + 1 

DO 42 K= 1 1 P 1 , N 

42 S = S + A { II , K ) * X ( K , J ) 

43 X(H»J5 = (X(II,J)-S)/A<II,II) 

KER= 1 

GO TO 15 
50 KER=2 

‘75 CONTINUE 
END 
END 

430 



APPENDIX II 



GRAPHS 

OP SIMULATED SYSTEM TRAJECTORIES 
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